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Preface 


This is the final report on this project which was concerned with the analysis of 
cavity-backed antennas and more specifically spiral antennas. The project was a continua- 
tion of a previous analysis, which employed rectangular brick elements, and was, thus, re- 
stricted to planar rectangular patch antennas. A total of five reports were submitted under 
this project and we expect that at least four journal papers will result from the research de- 
scribed in these reports. The abstracts of the four previous repons follow this page. The 
first of the reports (02891 8- 1-T) is over 75 pages and describes the general formulation 
using tetrahedral elements and the computer program. Report 02891 8-2-T was wntten af- 
ter the completion of the computer program and reviews the capability of the analysis and 
associated software for planar circular rectangular patches and for a rectangular planar spi- 
ral. Measurements were also done at the University of Michigan and at Mission Research 
Corp. for the purpose of validating the software. We are pleased to acknowledge a partial 
support from Mission Research Corp. in carrying out the work described in this report. 
The third report (028918-3-T) describes the formulation and partial validation (using 2D 
data) for patch antennas on a circular platform. The 3D validation and development of the 
formulation for patch antennas on circular platforms is still in progress. 

The fourth report (02891 8-4-T) is basically an invited journal paper which will ap- 
pear in the J. Electromagnetic Waves and Applications in early 1994. It describes the ap- 
plication of the finite element method in electromagnetics and is primarily based on our 
work here at U-M. 

This final report describes the culmination of our efforts in characterizing complex 
cavity-backed antennas on planar platforms. The report describes for the first time the 
analysis of non-planar spirals and non-rectangular slot antennas as well as traditional planar 
patch antennas. The comparisons between measurements and calculat ions are truly im- 
pressive. Another unique aspect of this work is the incorporation of the FFT as part of the 
BiCG solver by overlaying a structured triangular mesh over the unstructured mesh. The 
implementation of this BiCG-FFT solution algorithm is important in minimizing the CPU 
and storage requirements. This final report will be submitted for publication in a refereed 
journal. 
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A FINITE ELEMENT BOUNDARY INTEGRAL 
FORMULATION FOR RADIATION AND SCATTERING 
BY CAVITY ANTENNAS USING TETRAHEDRAL 

ELEMENTS 


J. Gong , J.L. Volakis, A. Chatter jee and J.M. Jin 
Radiation Laboratory 

Department of Electrical Engineering and Computer Science 
The University of Michigan 
Ann Arbor, Michigan 4 S 109-21 22 


ABSTRACT 

A hvbrid finite element-boundary integral formulation is developed us- 
ing tetrahedral/triangular elements for discretizing the cavity /aperture of 
microstrip antennas/arrays. The tetrahedral elements with edge based Un- 
ear expansion functions are chosen for modeUng the volume region and tri- 
angular elements are used for discretizing the aperture. The edge-based 
expansion functions are divergenceless thus removing the requirement to 
introduce a penalty term and the tetrahedral elements permit greater ge- 
ometrical adaptabihty than the rectangular bricks. The underlying theory 
and resulting expressions are discussed in detail together with some numer- 
ical scattering examples for comparison and demonstration. 
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Electromagnetic scattering and radiation from 
microstrip patch antennas and spirals residing 

in a cavity 

J. L. Volakis, J. Gong, and A. Alexanian 


Abstract 

A new hybrid method is presented for the analysis of the scat- 
tering and radiation by conformal antennas and arrays comprised of 
circular or rectangular elements. In addition, calculations for cavity- 
backed spiral antennas are given. The method employs a finite ele- 
ment formulation within the cavity and the boundary integral (ex- 
act boundary condition) for terminating the mesh. By virtue of the 
finite element discretization, the method has no restrictions on the 
geometry and composition of the cavity or its termination. Further- 
more, because of the convolutional nature of the boundary integral 
and the inherent sparseness of the finite element matrix, the storage 
requirement is kept very low at O(n). These unique features of the 
method have already been exploited in other scattering applications 
and have permitted the analysis of large-size structures with remark- 
able efficiency. In this report, we describe the method’s formulation 
and implementation for circular and rectangular patch antennas in 
different superstate and substrate configurations which may also in- 
clude the presence of lumped loads and resistive sheets/cards. Also, 
various modelling approaches are investigated and implemented for 
characterizing a variety of feed structures to permit the computation 
of the input impedance and radiation pattern. Many computational 
examples for rectangular and circular patch configurations are presen- 
ted which demonstrate the method’s versatility, modeling capability 
and accuracy. 
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A Finite Element-Boundary Integral Method 
for Cavities in a Circular Cylinder 

Leo C. Kempel and John L. Volakis 
December 22, 1992 


Abstract 

Conformal antenna arrays offer many cost and weight advantages 
over conventional antenna systems. However, due to a lack of rigorous 
mathematical models for conformal antenna arrays, antenna design- 
ers result to measurements and planar antenna concepts for designing 
non-planar conformal antennas. Recently, we have found the finite 
element-boundary integral method to be very successful in modeling 
large planar arrays of arbitrary composition in a metallic plane. Here- 
with, we shall extend this formulation to conformal arrays on large 
metallic cylinders. In this report, we develop the mathematical for- 
mulation. In particular we discuss the shape functions, the resulting 
finite elements and the boundary integral equations, and the solution 
of the conformal finite element-boundary integral system. Some valid- 
ation results are presented and we further show how this formulation 
can be applied with minimal computational and memory resources. 
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A class of hybrid finite element methods 
for electromagnetics: a review 

John L. Yolakis, A. Ckatterjee and J. Gong 
Radiation Laboratory 

Department of Electrical Engineering and Computer Science 
The University of Michigan 
Ann Arbor MI 48109-2122 


Abstract 

Integra] equation methods have generally been the workhorse for antenna and scat- 
tering computations. In the case of antennas, they continue to be the prominent 
computational approach, but for scattering applications the requirement for large- 
scale computations has turned researchers’ attention to near neighbor methods such 
as the finite element method, which has low O(N) storage requirements and is read- 
ily adaptable in modeling complex geometrical features and material inhomogeneities. 
In this paper, we review three hybrid finite element methods for simulating compos- 
ite scatterers, conformal microstrip antennas and finite periodic arrays. Specifically, 
we discuss the finite element method and its application to electromagnetic problems 
when combined with the boundary integral, absorbing boundary conditions and arti- 
ficial absorbers for terminating the mesh. Particular attention is given to large-scale 
simulations, methods and solvers for achieving low memory requirements and code 
performance on parallel computing architectures. 
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A hybrid finite element-boundary integral 
method for the analysis of cavity-backed 
antennas of arbitrary shape 


Jian Gong, John L. Volakis, A.C. Woo and H.T.G. Wang 


Radiation Laboratory 
Department of Electrical Engineering 
and Computer Science 
The University of Michigan 
Ann Arbor MI 48109-2122 


Abstract 

An edge-based hybrid finite element-boundary integral (FE-BI) 
formulation using tetrahedral elements is described for scattering and 
radiation analysis of arbitrarily shaped cavity-backed patch anten- 
nas. By virtue of the finite element method (FEM), the cavity ir- 
regularities, the dielectric super/substrate inhomogeneities and the 
diverse excitation schemes inside the cavity may be readily modeled 
when tetrahedral elements are used to discretize the cavity. On the 
aperture, the volume mesh reduces to a triangular grid allowing the 
modeling of nonrectangular patches. Without special handling of the 
boundary integral system, this formulation is typically applicable to 
cavity-backed antenna systems with moderate aperture size. To retain 
an O(N) memory requirement, storage of the full matrix due to the 
boundary integral equation is avoided by resorting to a right triangu- 
lar aperture grid on taking advantage of the integral’s convolutional 
property. If necessary, this is achieved by overlaying a structured tri- 
angular grid on the unstructured triangular grid and relating the edge 
field coefficients between the two grids via a narrow banded, sparse 
transformation matrix. The combined linear system of equations is 
solved via the biconjugate gradient (BiCG) method, and the FFT al- 
gorithm is incorporated in the solver to compute the matrix- vector 
product efficiently, with minimal storage requirements. 
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1 Introduction 

Microstrip antennas have been extensively investigated experimentally, ana- 
lytically and numerically for decades. By and large, numerical methods have 
been serving the engineers and researchers in the analysis and design of these 
conformal antennas for many years. Among them the moment method in 
conjunction with various integral equation (IE) formulations placed a major 
role [1-3]. However, IE methods are associated with field representations in 
which the appropriate Green's function for the specific geometry must be 
employed and this limits their versatility. Moreover, IE techniques are usu- 
ally formulated on the assumption of an infinite substrate, a model which 
obviously deviates from the practical configuration leading to inaccuracies 
for larger bandwidth antennas. Furthermore, in the context of IE methods, 
antenna excitations are represented using simplified models that differ more 
or less from the actual configurations. Also, due to the singularity of the cur- 
rent distribution near the patch-probe junction(s), special measures must be 
taken [4], not to mention additional IE complexities due to possible substrate 
anisotropies or inhomogeneities in the antenna substructure. In contrast, the 
hybrid Finite Element-Boundary Integral (FE-BI) technique alleviates these 
difficulties and this was demonstrated recently when the method was applied 
to rectangular patch antennas [5]. 

In this paper, we present an edge-based hybrid finite element-boundary 
integral formulation using tetrahedral elements for a characterization of ar- 
bitrarily shaped cavity-backed antennas. An example of such a configuration 
is shown in Figure 1, where a cavity is recessed in a metallic ground plane 
enclosing the FEM volume and the antenna elements on the aperture may 
be excited by different schemes, such as a simple probe, a magnetic frill gen- 
erator, a practical coaxial cable, microstrip lines, slots or a CPW line. In the 
context of the FEM, the cavity is first discretized into a number of tetrahed- 
ral elements that naturally reduce to triangles on the cavity’s aperture. For 
nonrectangular patches this triangular gridding is, in general, non-uniform 
and the exact boundary integral formulation based upon this mesh applies 
to any patch shape. As a result, the hybrid FE-BI technique is capable of 
modeling arbitrarily shaped cavity-backed antenna configurations, different 
substrate inhomogeneities, anisotropies, as well as various practical excita- 
tion schemes. 

As is well known, the boundary integral (BI) equation subsystem leads 
to a fully populated matrix whose size is determined by the number of aper- 
ture mesh edges. For large apertures, this analysis becomes impractical in 
terms of storage and computation time requirements, and to overcome this 
inefficiency, a uniform zoning of the aperture is required as shown in Fig- 
ure 2. By resorting to the structured mesh, the boundary integral matrix 
can be cast into a discrete convolutional form, thus permitting the computa- 
tion of the matrix- vector products via the discrete Fourier transform (DFT) 


2 


avoiding a need to store the full BI matrix. This memory saving scheme has 
already been applied to IE solutions involving rectangular surface grids [5,6] 
and in this paper, we describe how the BiCG-FFT solution is implemented 
for triangular meshes. The differences between the rectangular and triangu- 
lar meshes are also described and results are presented which demonstrate 
the method's vesatility in computing the scattering and input impedance of 
various non-rectangular printed antennas. 


2 Formulation 

In this section, we present the edge-based hybrid FE-BI formulation using 
variational principles, where the matrix algebra notation is employed so that 
we can readily extend the formulae to the general anisotropic case. We begin 
with the vector wave equation, 

V x x E j - k*e r E = -jk 0 Z 0 J' - V x (1) 

where J, and M, represent interior electric and magnetic current sources 
within the cavity V, e r and /z r denote, respectively, the relative permittivity 
and permeability, and k 0 is the free space wave number. The solution of (1) 
under certain boundary conditions is equivalent to extremizing the functional 

F(E,H) = ^ JJJ V {( v x E) • “(V x E) — &o £ rE • E j dv 
+ jk 0 Z 0 JJ s E xH-dS 

+ JJJ v E • (jk Q Z 0 J, + V x dv (2) 

where the surface S emcompases the cavity aperture excluding the portion 
occupied by the antenna elements. To proceed with the extremization of F, 
it is necessary to write H in terms of E. This can be accomplished through 
the surface equivalence principle. Following this procedure, the magnetic 
field H in the region exterior to the cavity may be expressed in terms of the 
equivalent magnetic surface current 

M = 2E x z (3) 

whose support is only over the extent of the surface S and the factor of 2 is 
due to image theory. Specifically, we have 

H(r) = H' + H r - 2 jk 0 Yo jj ( I + ^ VV ) G «(r, r ') ' ( E x *) dS ‘ ( 4 ) 
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where H‘ is the incident field, if any, from the exterior region, H r is the 
reflected field due to H' in the absence of the cavity, I is the unit dyad, 
and Go(T.r') is the free space Green’s function with r and r 1 denoting the 
observation and integration points. Substituting (4) into (2) yields 

F( E) = ;|j(V X E).i(Vx E )-^E.E}^ 

+ 2jk 0 Z 0 JJ s (ExH i )-zdS 

+ jjj^ E • (jfc 0 Z 0 J, + V x j-M t - J dv (5) 

- 2 kl jj s ( E x £) • | JJ s ( E x z) • (l + ^ w) Go(ry) dS 1 } dS 

which is now only in terms of electric field. In deriving this, we used the 
relation n x H r = h x H‘ which is valid when H r is measured on the ground 
plane. 


2.1 Cavity Volume Modeling 

In proceeding with the discretization of (5), it is convenient to re-express it 
as 

F = Fv + Fs (6) 


where Fv denotes the volume integral contributions and similarly Fs accounts 
for the surface integral contributions. The cavity volume is subdivided into 
N tetrahedral elements V e (e=l,2,...N), and within each tetrahedron the field 
is expanded as 

E = [V)l{E} e (7) 

with 


Me 

{Vu} 


{E}e 


[{Vr}{V y }{V z }) e 


Kl \ 

V u2 

V K6 J 


u = x,y,z 


(l) 


\Ee/' 


( 8 ) 


in which V ut is the u (u = x,y or z) component of the volume vector basis 
functions along the ith edge. The unknown vector {E} e has six entries, one 
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for each tetrahedron edge. (In this paper, we use square brackets for matrices 
and curly brackets for vectors). Inserting (7) into (5), and taking the first 
variation of Fy with respect to {£} e , yields 

6F V =L{[- 4 M £ }' + < A 'M (9) 


w 


here 


M], 

{!<}. 


= jJJ v IFdV],\DV]J - *o£,[l'WVirj dv 




[V)A 


jk 0 Zo 



\dv ir = 


TyW-T: K) 

IW-Ty™ 



( 12 ) 


To carry out the above integrations, it remains to introduce the volume 
expansion or shape functions V e . For our implementation we employed the 
linear edge- based shape functions for tetrahedral elements given in [7,8]. 


2.2 Aperture Modeling 

To discretize the surface integrals in (5), the aperture is subdivided into 
triangular elements since these correspond to the faces of the tetrahedrals. 
Within each triangle, the field is represented as 


e = {s)Ue.}, 

where 

[si. = HS.HS,}], 

( s " \ 

{•?u} = I $u 2 I i u = x, y 

\ Su3 J 

( E. x \ 

{E 3 } t = E s2 

\E*). 


(13) 


( 14 ) 
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where S ut is the u{u = x.y) component of the surface vector basis functions 
along the ?th edge. On substituting (13) into the surface integrals in (5) and 
taking the first variation of Fs with respect to {Ej},.. we obtain 

= + < 15 ) 


where 

[B} e = [I ff |-2*g[Se][; 

JJS' JJs 1 ' I 


+ 






•G 0 (r, r') dS dS 1 


(16) 


and 


{L} e =j2k 0 Z 0 JJ s [S e ] 



dS 


(17) 


Note that in (16) the elements of the array [S e ] are functions of the observa- 
tion vector r, whereas the elements of [5 e ]^ are with respect to the integration 
point r'. A suitable set of linear edge-based surface basis functions is 


S,(r) 


- l -j-z x (r - r,)e(r) r € S e 
0 otherwise 


(18) 


In this expression, /, denotes the length of the ith edge and is the position 
vector of the vertex opposite to the ith edge. Since each edge shares two 
triangles, one is defined as the plus and the other as the minus triangle. 
Therefore, e(r) is given by 


t(r) 


1 r € S e + 
-1 res e - 


(19) 


where S e = S+ + S~. The constant A e in (18) denotes the area of the plus 
or minus triangle depending on whether r € S+ or r € S~ . We note that 
S;(r) x z yields the basis functions used by Rao, etc. [9] in their moment 
method solution of boundary integral equations. 


2.3 System Assembly 

To construct the final system for the solution of the electric field components 
we combine (9) and (15), and after assembly we obtain the system 

{[*]{£} + f/0) + {[£]{£.} + {ill = o ( 20 ) 

In this, {A'} and { L } are the excitation vectors due to the interior current 
sources and the exterior excitation, respectively. The unknown electric field 
vector {E} consists of all field expansion coefficients with respect to the 
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element edges except those coinciding with perfectly electrically conducting 
(PEC) walls, PEC antenna element(s) or PEC pins inside the cavity. Finally, 
the vector {F,} represents the unknown surface fields whose entries are part 
of those in {£} with their corresponding edges on the aperture. Explicit 
expressions for the matrices and vectors in (20) are 


u 



A 7 

L\ 


III (V x V?) • (V X V') - kle r V' ■ V<) dv (21) 

2*o JJ S [i x SJ(r)] {jl [Sf(r') x z] C 0 (r,r')ds'j ds (22) 
+ 2 JJ $ ' ( S j( r ) x 5 ) | ff s V 5 ' [ S f( r ') x f ] G , 0 (r,r')ds , | ds 

IIi V r{j W + V,(i M ')(* (23) 

2jkoZo IIs s r [ H ' x "] ds (24) 


when specialized to inhomogeneous, isotropic dielectric cavity fillings. It is 
evident, from (21) and (22), that [A] and [B] are symmetric as a result of 
the assumed isotropic medium and reciprocity. In addition, [A] exhibits high 
sparsity due to the FEM formulation whereas [B] is fully populated. Two 
approaches may be followed in carrying out the solution of the combined sub- 
systems when an iterative solver is employed such as the biconjugate gradient 
(BiCG) method [10]. These two approaches differ in the manner used for the 
evaluation of matrix-vector products called for in the BiCG algorithm. One 
could sum the coefficient matrices [A] and [B] by adding up the correspond- 
ing matrix entries prior to the execution of the BiCG algorithm, or instead 
the resulting vectors may be summed up after carrying out the matrix- vector 
products. We observed that the first approach is more efficient in terms of 
computation time after reordering the combined matrix and storing only the 
non-zero elements. This is because in the context of this scheme the combin- 
ation of the two matrices is performed only once outside the BiCG iteration. 
However, the second approach is compatible with the BiCG-FFT scheme, 
when the FFT algorithm is employed to exploit the convolutional property 
of the integral operator, thus eliminating a need to explicitly store the en- 
tire BI matrix. Below, we discuss the implementation of the matrix-vector 
product of the boundary integral system for the BiCG-FFT solution. 


2.4 Implementation of the Boundary Integral Matrix 
Vector Products Using FFT 

To describe the execution of the matrix-vector product in the context of the 
BiCG-FFT solution scheme we refer to Figure 2. In this figure we display the 
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uniform triangular grid which either overlays the irregular aperture grid or 
is simply the aperture grid which the volume tetrahedral mesh is by design 
reduced to during the discretization process. We recognize that the proposed 
triangular grid consists of equal right triangles and thus it involves three 
different classes of edges (class 1, 2 and 3). These include the x-directed, 
(/-directed and the diagonal edges, all of which are uniformly spaced. For 
the FFT implementation each class of edges is independently numbered in 
accordance with their geometric location. Specifically, the ith class will carry 
the numbering (m, n) if the edge is the mth along the x direction and the 
nth along the y direction. The indices (m,n) take the values 

m = 0, 1, 2, ... , M' 
n = 0,1,2,..., N' 


with i — 1 for the y-directed edges, i = 2 for the diagonal edges and i = 3 
for the x-directed edges. Consequently, we find that 


M' 


' At - 2 i = l 
< M - 1 i = 2 
M- 1 i = 3 


N l = l 


{ N- 1 
N- 1 
N - 2 



(25) 


where M and N denote the numbers of elements along the x and y directions, 
respectively. 

To perform the integrations for the evaluation of the boundary integral 
matrix elements, it is now convenient to rewrite the basis functions (18) in 
terms of the new indices (m,n). We readily find that the edge-based basis 
functions associated with each of the aforementioned class of edges can be 
rewritten as 


SL(^y) 


Sl n (x,y) 


s 3 mn (^y) 


l 

Ax 


(nAy - y)x + (x — mAx)y (x,y)eS* 

{y- (n+ l)Ay)x + ((m + 2)Ax - x)y (x,y)eS; (26) 
0 otherwise 


VjAxg + (Ay) 2 f 
AxAy 1 


(nAy - y)x -(- (x - (m + l)Ax)y 
(y - (n -(- l)Ay)x -|- (mAx - x)y 
0 


f- 9 + 

es; (27) 
otherwise 



((n + 2)A y - y)x + (x - (m + l)Ax)y 
(y - nAy)x + (mAx - x)y 
0 


(x, y) G S* 
(x,y) € S~ 
otherwise 


(28) 


w r here the superscripts refer to the edge class. Each entry of the boundary 
matrix- vector product can now be calculated as 


3 

{BI subsystem) = [B] {E 3 } = 

i=i 


m'=0 n^O 


(29) 


in which (m,n) are the geometric location indices for the ith class edges of 
observation elements whereas ( m',n ') are the same for the _/th class edges of 
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testing elements. Thus, the specification of the indices i , m and n completely 
defines the entry k, = nM ’ + M of the column resulting after the execution 
of the boundary matrix- vector product. It is readily found that 


JJ. IL s '"" ■ s -'"' Go(r ’ r>) iz dy di ' dy ' 


(30) 


+ 


(AxAy ) 2 


jj Jj s e,{r) e } {r) l,l : G 0 {r,r') dx dy dx' dy' 


with 





^/(Ax ) 2 + (Ay ) 2 

Ax 



(31) 


More importantly, it can be shown that the BI subsystem exhibits the con- 
volutional property B'^ m , n , = £ ( t_ m ., n _ n .) and thus we can rewrite C 29 ) 
as 3 

[B}{E,} = j2 BiJ * £j ( 32 ) 


where the * denotes convolution. It is now seen that the computation of 
the boundary matrix- vector product can be performed by employing the 2-D 
discrete Fourier transform (DFT) thus avoiding a need to store the BI matrix 
other than those entries which are unique. When the symmetry property of 
£(m-m',n-n') 5s also invoked, implying 




i — n ' 1 


— R -? 1 

^ (m* —m^n 1 —n) 


(33) 


it is concluded that the total non-redundant entries in the BI matrix are 


3 3 

iv P = EEW + M) - 1 ) 


«=i j=i 


(34) 


This should be compared to the (£?_, M'N'J entries whose storage would 
normally be required if the BI system was not cast in convolutional form. We 
remark that N p is nevertheless equal to twice the number of entries required 
for uniform rectangular grids [6] for one class of edges. To avoid aliasing, it 
is necessary that B^ m _ m , n _ n ,^ = B^(m,n) be cast in a 2-D array which has 
the usual periodic form and zero padding may also be required to make use 
of the standard FFT routines. Specifically, the matrix-vector product (29) is 
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executed by using the MFTxNFT array 


~ ■"> 0 < m < M' 

o < n < s i 


B' } (-m, -n), 


MFT - M i + 1 < m < MFT 
0 < n < -V J 


B' p 3 ( m,n) — <J (fh n — l - NTT) 0 < m < M 

b {m,n i in r l ), NFT _ + \ < n < NFT 

^ MFT — M' + 1 < fh < MFT 

B' 1 (m - \ - MFT,n - 1 - NFT), j^ FT _ + i < n < NFT 


0 otherwise 

with the corresponding field vector given by 

_ ?/ _ f E J (m,n) 0 < m < M\ 0 < n < jV j 

Cp(m,n) = | o otherwise 


(35) 

(36) 


and MFT and NFT must be powers of 2 if a radix 2 FFT is used. 

In the BiCG-FFT algorithm the BI subsystem vector is computed as 

3 

{BI subsystem} = jDFT -1 |DFT{Bp J } ■DFT{£'p}|| (37) 

.=1 


The presence of the operator 5 indicates the necessary reordering of the 2-D 
array which results after the inverse FFT operation into a single column with 
the proper indexing for addition to the FEM subsystem. 


3 Mesh Overlay Scheme for Nonrectangular 
Patches 

As described above, the BiCG-FFT solver requires uniform aperture gridding 
so that the BI subsystem can be put in block circulant form. This can be 
always achieved during mesh generation whenever the patches are rectangular 
in shape or in case of radiators which are placed at some distance (usually 
small) below the aperture. However, for circular, triangular, or other non- 
rectangular patches on the aperture, it is not possible to construct a uniform 
mesh using the mesh generator. Typically, the aperture mesh is necessary 
to comform to the patch shape, leading to an unstructured surface grid. In 
this case, to make use of the efficient, low memory BiCG-FFT algorithm, an 
approach is to overlay on the unstructured aperture grid another coincident 
structured grid as shown in Figure 3. The boundary integral subsystem is 
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then constructed by using the overlaid uniform grid whose edge fields can be 
related to those on the unstructured grid via a sparse transformation matrix. 
That is, it is necessary to append to the system (20) the relations 


{£.}„ = [TF]{E,) nu 

{£>}„„ = |r B ]{£.), (38) 

where the subscripts u and nu refer to the field coefficients of the uniform 
and non-uniform aperture grids, respectively. Also, [7>] and [Tb] refer to the 
forward and backward transformation matrices, respectively, with A’ u and 
N nu denoting the numbers of the uniform and non-uniform mesh edges on 
the cavity aperture. 

To derive the elements of [T>], we begin with the expansion (13) and 
enforce it at three points on each edge belonging to the uniform grid. We 
conveniently place these three points at the center and ends of the edge (see 
Figure 4). Given the fields at these points, we can approximate the field 
along the (m, n) edge of the uniform grid using the weighted average 

(E u ) (m ,„) = -e u ■ l — - 53 E* u (r en di) 

1 ^mid 

+ -T7— £ En.!^) 

^ ' mid fc = \ 

+ T, 5— £ E‘,(r.„d,) (39) 

^end 2 k = 1 J 


in which e u denotes the unit vector along x,y or the diagonal, depending 
on the class of edge being considered. The quantities E* u represent the 
fields in the non-uniform grid triangles with the superscript k being a sum 
variable in case Tend,, r en d 2 or Tmia specify a point shared by more than one 
triangle. Obviously, A^, A^mid and N end3 denote the number of non-uniform 
grid triangles sharing the node at r en d,, r m id and r e „d 2 , respectively, and will 
typically be equal to unity. 

After assembling (39) into (38) we find that the elements of the forward 
transformation matrix are given by 


i | i ^ end i 3 

(Tf)ij = - 53 13 c O/S < (r e „d 1 ) 

^mid 3 


1 ‘ ’mid 'J 

+ TT — 13 e «j<Sr(rmid) 

■‘’mid k=l (= i 
■j ^endj 3 

T rrr 53 X3 C 0<Sf(r e nd 2 ) 

2A/end 2 k=l t=l 


(40) 
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in which 


e . = 1 •>=•>'. 

Jt 1 0 otherwise 

and the global indices i and j correspond to the ith uniform grid edge and 
the j th non-uniform grid edge. The subscripts j ( is the global index used 
in numbering the non-uniform grid edges, whereas the subscript ( (= 1, 2 
or 3) is the local edge index used in the definition of the basis functions 
(see (13)). We remark that the explicit computation of the transformation 
matrix elements results in a substantial programming simplification because 
it avoids the usual assembly process. Following the same procedure, we can 
obtain the elements of the backward transformation matrix. Assuming each 
uniform grid edge traverses on three or less non-uniform grid triangles, the 
non-zero entries in each row of [TV] will be 9 or less. However, they can 
reach a maximum of 18 if the midpoint and endpoints reside on an edge of 
the non-uniform grid. The maximum non-zero entries in each row of [7b] 
will be 15, but the typical number will be much less. 


4 Numerical Considerations 

Based on the presented FE-BI formulation, a computer program was written 
for the analysis of the radiation and scattering by cavity-backed patch anten- 
nas of arbitrary shape. The antenna geometry is supplied to this program in 
an input file which as a minimum, it must contain lists of (a) the nodes and 
their (x,y,z) coordinates, (b) the nodes forming each tetrahedron, (c) the 
nodes on the cavity aperture, and (d) the nodes on metallic boundaries. 
For arbitrary antenna geometries, it is necessary to employ a sophisticated 
volume mesh generation package and a number of these are available com- 
mercially. Typically each of these packages generates a “Universal file” which 
can be readily preprocessed to extract the aforementioned input lists. 

A major effort was devoted in writing the program in a manner which 
minimizes the storage and computational requirements. Specifically, the BCs 
on the metallic surfaces are enforced a priori to obtain a system which in- 
volves only nonzero field components. The sparse finite element matrix was 
stored as a single array of length N„N m , where N v is the total number of 
unknowns within the cavity volume and iV nz denotes the maximum number 
of nonzero row entries. The BI matrix was stored in different ways depending 
on whether the FFT was to be employed for the evaluation of the matrix- 
vector products. If the BiCG solution was to be carried out without the 
FFT, then the N, x N s BI integral matrix was added to the FE array res- 
ulting in a 1-D array about N v iV nz + Nj long. For slot antennas, including 
cavity-backed spirals, and moderately sized systems, it was found preferable 
not to use the FFT, thus avoiding any discretization errors. In that case 
the generation of a single combined FE-BI matrix before execution of the 
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BiCG algorithm reduces the computational requirements because a number 
of operations associated with the repeated combinations of the FE and BI 
subsystems within the BiCG iteration is avoided. 

When the FFT is to be used as part of the BiCG solver the FE and BI 
matrices must be kept in separate arrays throughout the execution process. 
In this case the FE matrix is again stored as a single array and similarly 
the non-redundant elements of the BI matrix are stored in a single array 
of length 9.Y(2 M - 1). The factor of 9 is due to the three classes of edges 
and as usual M and N denote the number of elements along the i and y 
directions, respectively. Because of the storage and computational efficiency 
of BiCG-FFT algorithm, it is necessary to resort to uniform aperture grids 
for conformal antennas involving substantial number of aperture edges. Of 
course, one should always use uniform triangular grids when the patches 
on the aperture are rectangular or if the array is simply covered by a su- 
perstate. In the case of circular patches it will be necessary to overlay a 
structured triangular grid over the unstructured grid generated by the mesh 
generator. This must be done in the preprocessing stage and should be taken 
into account when constructing the FE matrix. For scattering computations, 
the overlay of the structured grid is almost always the preferred approach 
because it does not, generally, compromise the accuracy of the computed 
scattering cross section. However, for antenna parameter computations, the 
interpolation scheme between the structured and unstructured grid edges 
may be of concern, depending on the specific antenna geometry. Generally, 
circular slots and spirals should be treated without resorting to structured 
grids and to our experience this does not cause large computational burden 
because these antennas are associated with small apertures. In the case of 
circular patches, the structured grid was not seen to compromise the com- 
putational accuracy. Of course, conclusions based on one type of antenna do 
not necessarily apply to others, and thus the suggested alternatives must be 
examined separately for each antenna before choosing one approach over the 
other. Of importance here is that the formulation is suitable for modeling 
any antenna shape and feed structure. 


5 Results 

We present below some representative numerical results for the purpose of 
validating and demonstrating the robustness of the tetrahedral formulation 
for scattering and radiation by different configurations of cavity-backed an- 
tennas. In each case the computed results via the FE-BI method are com- 
pared with reference measured or calculated data. 

BiCG-FFT algorithm validation: To validate the BiCG-FFT implement- 
ation for the boundary integral matrix-vector products described in section 
2.4, we consider the scattering by a circular patch placed 0.325 cm be- 
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low the aperture of a rectangular cavity (recessed in a ground plane) of 
7.8 cm x 5.2 cm x 0.975 cm. The geometry is illustrated in Figure 5 and 
we remark that the patch was purposely submerged so that the structured 
grid using a mesh generation package was available on the aperture. A bi- 
static radar cross section (RCS) pattern is shown in Figure 5(a) as computed 
with and without use of the FFT in the BiCG algorithm for matrix-vector 
products. As seen, the two solution procedures give identical scattering pat- 
terns for two different cavity dielectric filling. 

Overlay scheme validation: To validate the structured/unstructured mesh 
overlay scheme discussed earlier, we calculate the scattering by a circular 
patch (0.1 cm radius) placed on the aperture of a rectangular cavity 0.4 cm 
per side and 0.1 cm deep. The bistatic radar cross section (RCS) patterns are 
shown in Figure 5(b) as computed with and without use of the overlay scheme 
for normal and 30° off ground plane incidence. The dielectric fillings were 
chosen as e T = 4 and e T = 2 for normal and oblique incidence, respectively. Of 
course, the FFT was only employed in connection with the overlaid structured 
mesh. 

Scattering by a circular patch: Figure 6 illustrates a circular patch resid- 
ing on the surface of a 0.406 cm thick substrate having a relative dielectric 
constant of e T = 2.9. The patch’s diameter is 2.6 cm and the substrate is 
enclosed in a circular cavity 6.292 cm wide. This cavity and the patch were 
recessed in a low cross section body for measuring its RCS. A comparison 
of the measured and calculated backscatter cree RCS as a function of the 
frequency is shown in Figure 7. For this computation the direction of the 
incident plane wave was 60° from normal, and as seen the agreement between 
measurements and calculations is very good throughout the 4-9 GHz band. 
Input impedance measurements and calculations for the same patch are dis- 
played in Figure 7. The feed in this case was placed 0.8 cm from the patch’s 
center and it is again seen that the measurements and calculations are in 
good agreement. 

Scattering by a square archimedean spiral: Figure 8 displays a square 
archimedean spiral. A square rather than a circular spiral was selected be- 
cause it allowed comparisons with calculations based on a different technique. 
The subject complimentary square spiral consists of strip arms, each of width 
0.09375 cm, and is placed (free-standing) on the aperture of a square air-filled 
cavity 2.812 cm per side and 0.9375 cm deep. The <f> = 45° plane <To$ and 
bistatic (incident plane wave was at 30° off normal) RCS patterns for 
this structure are shown in Figure 8. It is seen that the (Tee + RCS com- 
pares well with the corresponding data based on the finite difference-time 
domain (FD-TD) method [13]. No FD-TD data were available for the <r^ 
RCS pattern but as expected, the RCS return is much lower and vanishes 
at grazing. 

Annular slot impedance: Figure 9 shows a narrow circular (0.75 cm wide) 
annular slot situated in a circular cavity 24.7 cm wide and 3 cm deep. Be- 
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is very small for this application and as a result there is no need to invoke 
the FFT in the BiCG algorithm. The FE-BI method is basically quite ef- 
fective in modeling small aperture configurations without a need for special 
computational considerations. Input impedance calculations as a function 
of frequency for this radiator, excited by a probe placed across the slot, are 
shown in Figure 9, and agree well with the values calculated via a modal- 
boundary integral method [14]. For these calculations, the frequency was 
swept from 700-1000 MHz. The dielectric constant of the material filling the 
cavity was set to t T = 1.35 as in [14] and this is an effective value to account 
for the presence of a dielectric slot cover used as part of the measurement 
model for holding the plate. 

Radiation by a one-arm conical spiral: We considered the modeling of this 
radiator to demonstrate the geometrical versatility of the FE-BI method. A 
configuration of the spiral radiator and surface mesh is illustrated in Figure 10 
from different views. The top and bottom edges of the strip forming the 
spiral follow the lines p = 0.0503A exp[0.221(<^ ± 2.66)], z = a± exp(0.221<p), 
where (p,<f>.z) denote the standard cylindrical coordinates, a± are equal to 
0.0832A and 0.0257A, respectively, and 0 < <f> < 2k. This spiral arm resides 
on an inverted cone (9.24 cm tall) whose bottom cross section has a diameter 
of 1.68 cm and the top cross section has a diameter of 21.78 cm. For our 
calculations A = 30 cm (/ = 1 GHz) and the spiral was situated in a circular 
cavity 10.01 cm deep. The computed E# radiation pattern taken in the 
(t> = 90°-plane, using a probe feed at the cavity base, are given in Figure 11. 
It is seen that the E# principal plane pattern is in good agreement with the 
data given in [12]. However, the E$ pattern differs from the measured data 
primarily because of the circular cavity included in the analytical model. The 
latter was not part of the measurement configuration which consisted of the 
spiral antenna on a large circular plate. 

6 Conclusion 

We presented a hybrid finite element-boundary integral (FE-BI) formula- 
tion which incorporates linear tetrahedrals. The method was specifically 
developed for the radiation and scattering analysis of cavity-backed printed 
antennas, where the FEM is used for modeling the cavity region and the BI 
equation acts as a global boundary condition for terminating the mesh on 
the cavity aperture. The FE-BI formulation is particularly suited for the 
analysis of complex configurations and much emphasis was given here in de- 
veloping a solution technique requiring O(N) storage in spite of the resulting 
full BI subsystem. The latter was achieved by making use of the convolu- 
tional property resulting from the structured mesh, thus permitting use of 
the FFT in the BiCG solver for computing the matrix-vector products. For 
scattering calculations associated with large aperture structures, use of the 
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FFT proved essential in minimizing the computational requirements. 

A number of patches, slots as well as planar and non-planar spiral an- 
tennas were analyzed for the purpose of demonstrating the versatility and 
accuracy of the FE-BI technique. Certainly, the need to use a sophisticated 
mesh generation package is deterrent to the application of the technique for 
the analysis of simple antenna configurations. However, this is unavoidable 
when dealing with complex geometries and, moreover, the pervasive use of 
such commercial packages on desktop computers makes the technique quite 
attractive. 
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Figure 2. Illustration of the uniform (structured) right triangular grid. 
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Figure 3. Overlay of a structured aperture mesh over an unstructured mesh, 
shown here to conform to a circular patch. 
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Figure 4: Illustration of the parameters and geometry used in constructing the 

transformation matrix elements between the structured and unstructured mesh. 
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Figure 5(b). Comparison of the bistatic RCS patterns for a circular patch of radius 
0.1 cm (residing in a 0.4 x 0.4 x 0.1 rectangular cavity) as calculated 
with and without the mesh overlay scheme. The frequency was f =20 
GHz and the shown patterns correspond to two incidence angles, one 
at 6 = 0° (e r = 4) and the other at 6 = 60® (e r = 2) from normal. 
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Figure 6. Comparison of the computed and measured <r„ backseat ter RCS as 
a function of frequency for the shown circular patch. The incidence 
angle was 30® off the ground plane. 




Figure 7. Comparison of the computed and measured input impedance for the 
circular path shown in Figure 5. The feed was placed 0.8 cm from the 
center of the patch and the frequency was swept from 3 to 3.8 GHz. 
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Figure 10. Illustration of the configuration and mesh of the one-arm conical spiral 
used for the computation of Figure 11 (the mathematical deflations 
of the spiral edges are given in section 5). 
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Figure 11. Comparison of the calculated Eg and Ej, radiation patterns taken in 
the 4> = 90°-plane with data in reference [12] for the one- arm conical 
spiral shown in Figure 10. 
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